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We present the interatomic force constants and phonon dispersions of graphite and graphene from the LCBOPII empirical 
bond order potential. We find a good agreement with experimental results, particularly in comparison to other bond order 
potentials. From the flexural mode we determine the bending rigidity of graphene to be 0.69 eV at zero temperature. 
We discuss the large increase of this constant with temperature and argue that derivation of force constants from 
experimental values should take this feature into account. We examine also other graphitic systems, including multilayer 
graphene for which we show that the splitting of the flexural mode can provide a tool for characterization. 
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1. Introduction 

The phonon spectrum of a crystalline solid provides in- 
formation on several important physical properties like 
sound velocities, thermal conductivity, heat capacity and 
thermal expansion. The phonon spectrum of graphite has 
been intensively studied experimentally [IHS] and theoret- 
ically [THTU] in the past and some models have also been 
shown to be accurate for the description of fullerenes [5] 
and of graphite slabs [10!. In the more recent past, Ra- 
man spectroscopy has proven to be of crucial importance 
also for the characterization of graphene and nanotubes 
as well as for graphitic nanostructures of lower symme- 
try, like bent tubes and graphene edges [TTHT5] . However, 
the many unusual structural aspects of graphene, like the 
observed ripples [H] , negative thermal expansion [T51 [T5] , 
edge reconstruction [17] and localized [18] and extended 
defects [HI [20] make it desirable to describe the energetics 
of carbon in different structural and bonding configura- 
tions beyond the harmonic approximation by means of a 
unique potential. Bond order potentials are a class of em- 
pirical interatomic potentials (EIPs) designed for this pur- 
pose [2TU23] . They aim at describing not only the struc- 
ture around equilibrium but also anharmonic effects [24] 
and the possible breaking and formation of bonds in struc- 
tural phase transitions like the graphite to diamond tran- 
sition where the character of the bonding changes from 
sp 2 to sp 3 [251 ES]. In view of this larger and exacting 
scope it may be expected that the phonon spectra derived 
from these potentials are not as accurate as those derived 
from models meant to describe a single specific situation. 
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However they allow to study, without further adjustment 
of parameters, all carbon structures, including the effect 
of defects, edges and other structural changes, also as a 
function of temperature. The purpose of this paper is 
to evaluate the force constants and phonon spectrum of 
graphitic structures derived from the Long-range Carbon 
Bond Order Potential (LCBOPII) [2U [27] and compare 
these results to experimental values, force constant mod- 
els and to the Tersoff [28] and Brenner [22] bond order 
potentials for carbon. 

The phonon dispersions of graphene and graphite have 
been measured experimentally [1-6] [29] , determined from 
ab initio calculations [5[ 15011521 and calculated from bond 
order EIPs [33J [M] . The ab initio results generally agree 
very well with experimental measurements whereas widely 
used EIPs such as the Tersoff [28] and Brenner [22] po- 
tentials give less accurate results [33] . particularly in the 
optical region. One reason for this is that the range of 
interatomic interactions in EIPs is limited for computa- 
tional efficiency whereas force constant models show that 
interatomic force constants (IFCs) up to fourth or even 
fifth nearest neighbours (NNs) must be included for ac- 
curate phonon dispersions [TJ [3T] . The second generation 
LCBOPII [37] EIP includes long-range interactions up to 
6 A which is well beyond fifth NNs in graphene and it is 
interesting to study their effect on the phonons in compar- 
ison to other approaches. 

In Section [2] we describe the computational method 
with emphasis on the anomalous dispersion of the flex- 
ural mode. In Section [3] we present the LCBOPII phonon 
dispersion of graphene and graphite, compare our force 
constants to other models and examine the role of specific 
force constants on the phonon dispersions. We devote Sec- 
tion [4] to the analysis of the bending rigidity and its tem- 
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Figure 1: Phonon frequency in cm . Left: graphene phonon dispersion from LCBOPII; Right: Graphite phonon dispersion from LCBOPII 
with experimental data, the inset is an enlargement of the low-frequency dispersion along the A — V line. The locations of the high symmetry 
points are M = 7r/\/3a(\/3, 1, 0), K = 47r/3a(l, 0, 0) and A = n/c(0, 0, 1) in the coordinate system defined in the inset of the left figure. The 
experimental data for graphite are from Ref. [2] (squares), Ref. \T\ (circles), Ref. [3] (triangles), Ref. [5] (diamonds), Ref. [4] (inverse triangles) 
and Ref. |29) (pentagons). 



perature dependence. In Section [5] we show the phonons 
of (10,10) nanotubes and show the relevance of low-energy 
phonons of multilayer graphene for their characterization. 

2. Methods 

The phonon dispersions are calculated by means of stan- 
dard lattice dynamics |35) . The interatomic force con- 
stants are calculated by evaluating, by central differences, 
the second derivatives (the IFCs) of the LCBOPII EIP 
with respect to atomic displacements around their equi- 
librium positions. The phonon frequencies at wavevector 
q, tu(q), are determined by diagonalizing the dynamical 
matrix 



1 



iq-R 



(i) 



R 



,k,k' , 



where <^'^ (R) is the force constant matrix, a, /3 being 
Cartesian indices, for two atoms k and k' in unit cells 
separated by a lattice vector R. In layered materials the 
lowest, out-of-plane, acoustic phonon branch (ZA) has a 
peculiar quadratic dispersion near the zone center with 
a coefficient determined by the bending rigidity k of the 
crystal. For graphite, it has been shown by Lifshitz [36 
that the dispersion has the following form: 



w Z a(<?) 



' C44 J 

P3D 



P3DC 



(2) 



where Pzd — 8mc/(\/3a 2 c) is the mass density, C44 is the 
shear elastic constant, c is the lattice parameter equal to 
twice the interlayer distance in ABAB stacked graphite, 



a is the in-plane lattice parameter and mc is the atomic 
mass of carbon. For graphene, the dispersion reduces to a 
purely quadratic form: 



wza(<?) = \ j 191 

PlD 



(3) 



where P2D = 4mc/(v / 3a 2 ) is now a two-dimensional mass 
density. 

3. Phonon dispersion 

Minimization of the LCBOPII cohesive energy with re- 
spect to the lattice parameters gives a = VSacc — 
2.4592 A for graphene and a = V3a cc = 2.4585 A, 
c = 6.7344 A for graphite. The phonon dispersions for 
graphene and ABAB graphite, calculated at these lattice 
parameters, are shown in Fig. [l] The branches are clas- 
sified as follows: L stands for longitudinal in-plane, T for 
transverse in-plane and Z for transverse out-of-plane polar- 
ization at the r point. An A refers to acoustic modes and 
an O to optical modes. The O' modes in graphite indicate 
out-of-phase oscillation of two equivalent atoms in neigh- 
bouring layers. The phonon dispersions of graphite and 
graphene are very similar due to the weakness of the inter- 
layer interactions compared to the strong covalent bonds 
binding the atoms in the layers. Consequently, most of the 
twelve branches in graphite are almost doubly degenerate 
with the exception of the out-of-plane branches below 400 
cm . 

Contrary to the two linear, in-plane acoustic LA and 
TA, modes the out-of-plane ZA mode has a quadratic dis- 
persion near T which is typical of layered crystals |37j . 
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Table 1: Graphcne phonon frequencies from LCBOPII at high symmetry points in cm . Experimental values for graphite are also listed: 
a Ref. g], 6 Ref. 0, c Ref. [2], d Ref. [29] (graphene), e Ref. [3], f Ref. [5] and " Ref. [I]. 



Mode 




r 




M 




K 




LCBOPII 


Experiment 


LCBOPII 


Experiment 


LCBOPII 


Experiment 


ZA 







265 


471 a , 465 b , 451 d , 485 s 


405 


482 d , 517 d , 530 e , 540 s 


TA 







713 


630 d , 631 s 


1033 


1010 s 


LA 







1282 


1290 c 


1153 


1194 c , 1224' 1 


ZO 


797 


861 b , 868 s 


540 


670 b , 631 s 


405 


588 d , 627 e , 540 s 


LO 


1563 


1590 b , 1582 d 


1290 


1323 c 


1153 


1194 c , 1224 d 


TO 


1563 


1590 b , 1575^, 1582 d 


1441 


1390 c , 1389 b 


1513 


1310 d , 1291 e 



The ZA mode is a bending mode, the two atoms in the 
unit cell move in phase in the z-dircction, which, at long 
wavelengths, bends the surface resulting in rippling of the 
graphene sheet. The softness of this mode also means that 
it plays a dominant role at low temperatures. Also the op- 
tical out-of-plane, ZO mode has a considerably lower en- 
ergy than the other optical branches due to the fact that 
atoms are much more free to move perpendicular to the 
plane than in the plane itself. At the -ftT-point, the TO/LO 
and the ZA/ZO modes are degenerate by symmetry. 

In Fig. [T]and Table[l]we compare the LCBOPII phonon 
spectrum to experimental results by high resolution elec- 
tron energy- loss spectroscopy (HREELS) [3] (5] |2"5] , in- 
elastic x-ray scattering [TJ [2] and inelastic neutron scat- 
tering [3]. The overall agreement with the experimental 
values is rather good, considering that the potential was 
not specifically fitted to reproduce the force constants of 
graphite. 

LCBOPII performs very well compared to the popular 
Tersoff [28] and Brenner [22] EIPs, for which the phonon 
dispersions were recently published [33 . The Tersoff EIP 
overestimates the LO and TO branches by nearly 40% and 
both potentials show large discrepancies with experiments 
in the in-plane acoustic branches which are very well re- 
produced by LCBOPII. The latter modes are of particular 
importance for the thermal conductivity in graphene |38] . 
The only deviation occurs at the M point for the TA 
branch, where LCBOPII overestimates the experimental 
value from Ref. pQ by 13%. The measurements from 
Ref. 5j show even higher frequencies for this mode but 
these may have been obtained from a sample of poor qual- 
ity as HREELS selection rules state that the TA mode 
should not be observable along the T — M line [1, 29 . Ab 
initio calculations also confirm the experimental results 
from Ref. pQ . From the slopes of the TA and LA modes 
we determine their respective sound velocities as 13.0 and 
20.7 km/s which compare well to the experimental values 
of 14.7 and 22.2 km/s [35]. 

The quadratic dispersion of the ZA mode is reproduced 
well by LCBOPII, but the frequency is underestimated. 
We argue that this mode is strongly temperature depen- 
dent as we discuss in Section [4j Also the ZO mode is 
found to be lower than all experiments, the difference be- 



ing about 8% at T. 

The low-energy dispersion of graphite for wavevectors 
parallel to the c-axis is shown in the inset of Fig. [I] In these 
modes, layers oscillate rigidly and frequencies are thus de- 
termined by the long-range interactions of LCBOPII. The 
two longitudinal, LA and LO', 'breathing' modes are in 
excellent agreement with experiments. This comes as no 
surprise since the compressibility of graphite was one of the 
parameters used in fitting the long-range interactions |27] . 
The lower, doubly degenerate, transverse branches, corre- 
sponding to the shearing motion of the layers are instead 
too soft compared to experiments as we discuss later in 
Section g] 

The description of the highest optical branches requires 
long-range IFCs up to fourth or fifth NNs [TJ [3TJ and are 
therefore the most difficult to reproduce with EIPs [33 . 
LCBOPII includes interactions in this range through its 
long-range potential, but these are merely pair interactions 
of Morse form which are too smooth to produce significant 
force constants. In particular, the flat dispersion of the TO 
mode along the M — K line differs from experiments and 
ab initio calculations. The difference reaches 15% at the 
.RT-point. Another missing feature is the overbending of 
the LO mode, namely the shift of the highest frequencies 
away from T. This overbending is believed to originate 
from strong electron-phonon interactions which lower the 
frequencies of the highest phonon modes at T and K [3] 
[40] . For these branches the Brenner EIP gives a behaviour 
similar to LCBOPII. 

To gain more insight in the origin of the discrepancies 
with experiments we compare in Table [2] the force con- 
stants of LCBOPII to IFC sets proposed in Refs. [UGHHM]. 
The coordinate system is chosen such that x is the co- 
ordinate along the line connecting two atoms, y is the 
in-plane coordinate perpendicular to this direction and z 
is the coordinate perpendicular to the plane. The IFCs 
between i-th NNs in these directions are respectively the 
bond stretching, <f$ , the transverse, 4>[l , and the out-of- 
plane, <pi , force constants. The general form of the force 
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Table 2: Force constants for graphcnc from LCBOPII compared to force constant models from Refs. 1. 31, 34 which were fitted to reproduce 
experimental results. 



i 

(Distance) 


Stretching, ^ J (eV/A 2 ) 
LCBOPII PQ [31] [34] 


In-plane, ^ (eV/A 2 ) 
LCBOPII [J [3T] [24] 


Out-of-plane, ^ (eV/A 2 ) 
LCBOPII | P i 


1 (1.42 A) 

2 (2.46 A) 

3 (2.84 A) 

4 (3.76 A) 

5 (4.26 A) 


26.60 25.88 25.58 25.57 
3.37 4.04 4.63 -2.55 
0.51 -3.02 -2.07 -2.07 
0.02 0.56 0.41 0.66 
0.00 1.03 0.00 0.00 


8.99 8.42 9.05 9.05 
-0.61 -3.04 -2.55 4.63 
-0.05 3.95 3.13 3.13 
0.00 0.13 0.34 0.31 
0.00 0.11 0.00 0.00 


4.73 6.18 6.17 6.17 
-0.75 -0.49 -0.51 -0.51 
-0.05 0.52 0.36 0.36 
0.00 -0.52 -0.32 -0.33 
0.00 0.17 0.00 0.00 



constant matrix for the i-th NNs in graphene is 

(4® o\ 

4? o , (4) 

V o o ^7 

where the off-diagonal, (f>^l, elements for i — 1,3,5 are 

(2) 

equal to zero due to the hexagonal symmetry, (f> y o( { — 1.48 

eV/A 2 and <j$ = O(10~ 6 ) eV/A 2 for LCBOPII. The IFCs 
from Rcf. [T] are obtained by fitting a fourth NN force con- 
stants model to experimental values, those from Ref. [34 
are derived from an extended Brenner EIP, also fitted to 
experimental values, and the IFCs from Ref. [31] are from 

(2) 

a fourth NN force constant model including a nonzero d 
of -0.57 eV/A 2 fitted to the dispersion obtained ab initio 
within the DFT-GGA approximation. 

From the comparison of Table [2] we see that the first 
NN IFCs are very similar to the fitted force constants of 
the reference models. However the decay of the LCBOPII 
force constants beyond first NNs is too rapid compared to 
the sets of IFCs which reproduce the experimental values 
accurately. 

To see how the in-plane phonon branches evolve if larger 
IFCs beyond first NNs are included we manually increase 
these force constants. We consecutively changed the in- 
plane IFCs from second to fifth NNs to match those from 
Ref. pp. Since <\>\ is already 23% lower than the fitted 
values we changed the out-of-plane IFCs from first NNs. 
They were matched to those obtained in Ref. [34] since 
their model includes interactions up to fourth NNs only 
which better resembles LCBOPII. 

The resulting phonon dispersions are shown in Fig. [2] 
The modification of the in-plane IFCs beyond first NNs 
clearly improves the phonon dispersion. Changes up to 
third NNs lower the frequencies of the transverse branches, 
particularly the TO branch along M — K and the TA 
branch at the TVf-point but perturbs the good agreement 
of the sound velocities of the linear modes. With in- 
plane IFCs changed up to fifth NNs the dispersion is in 
excellent agreement with experiments. This means that 
an improvement of only the long-range interactions of 
LCBOPII can considerably improve the accuracy of the 
potential for graphitic systems. However, this is not an 
easy task, since in the construction of LCBOPII short and 
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Figure 2: Left panels: modification of the in-plane force constants 
from the second to the indicated level of NNs to match those from 
Ref. pQ; Right panels: modification of the out-of-plane force con- 
stants from the first to the indicated level of NNs to match those 
from Ref. 134] . Notice the different scales. Red (solid) lines are the 
branches that are modified as a consequence of the change in force 
constants. The gray (dashed) lines are the original LCBOPII disper- 
sions. 



long range interactions are fitted simultaneously. For the 
out-of-plane branches the optical ZO branch is greatly im- 
proved by the increase of the first NN IFC but the im- 
portant quadratic behaviour of the ZA mode is lost. In- 
terestingly the quadratic dispersion is recovered only once 
fourth NNs interaction is included. 

4. Quadratic ZA dispersion and bending rigidity 

The bending rigidity k is the key quantity which charac- 
terizes the mechanical properties of membranes |41j . For a 
crystalline membrane like graphene it is intimately related 
to the quadratic ZA branch through Eq. ([3| . Reported val- 



ues of k vary from 0.79 to 2.13 eV [Ml El USS] . Besides 
the different techniques used to calculate n and the dif- 
ferent models of carbon interactions there might be other 
reasons for the confusing variety of reported values. First, 
when comparing the values of k for graphite and graphene, 
one should consider the bending rigidity per layer and not 
per unit cell since the latter results in a factor two differ- 
ence. In fact, since graphite has two graphene layers in the 
unit cell, the coefficients of the |<7| 2 term in Eqs. ^ and ^ 
differ by a factor two while (see Fig. [T]) the quadratic coef- 
ficients should be approximately equal. The second, more 
important reason, is that the bending rigidity of graphene 
has been found to be strongly temperature dependent in 
detailed Monte Carlo simulations [Ml H2]- Contrary to 
liquid membranes |41j . k increases with increasing tem- 
perature. This increase reaches roughly 40% already at 
room temperature [2H H2] ■ The temperature dependence 
of the bending rigidity implies that also the ZA phonon 
mode should depend on temperature which makes com- 
parison of zero temperature dispersion as presented here 
to the room temperature experimental values non straight- 
forward. 

From the fit of Eq. ([3| to the ZA dispersion along the 
r — M line we obtain for the bending rigidity n = (0.69 ± 
0.02) eV for graphene. The same procedure for graphite 
with Eq. ([2} yields k = (0.69 ± 0.01) eV per layer and 
C 44 = (5.8 ± 0.02) x 10 8 Pa. This value of C 44 is much 
lower than the experimental value of 5.03 x 10 9 Pa [3"9"j 
due to the too small corrugation energy of LCBOPII which 
gives a difference of about 1.5 meV/atom against about 10 
meV/atom [4"S1 ITT] bet ween AA and AB graphite stacking. 
As a consequence, the transverse modes of graphite for 
wavevectors parallel to the c-axis shown in the inset of 
Fig. [T] have too low frequencies. The bending rigidities of 
graphite and graphene per layer are the same, which is in 
agreement with the results from Ref. [32] who find almost 
equal bending rigidities per layer for graphene and bilayer 
graphene. 

The value of 0.69 eV is low compared to other stud- 
ies and also lower than the 0.82 eV reported earlier for 
LCBOPII at zero temperature |24j . This value was ob- 
tained by evaluating the elastic energy per unit area £ of 
carbon nanotubes. This energy is equal to £ = \/2nH 2 
where H is the curvature of the nanotube. The apparent 
discrepancy with the result from the phonons is due to the 
fact that in forming a nanotube from graphene both elastic 
and torsional energy occur and it is only in the limit of very 
large nanotubes that the torsion energy can be neglected. 
We calculated the elastic energy, and the corresponding 
bending rigidity, for several nanotubes with radii from 11 
to 70 A with and without the inclusion of the torsion term. 
The results in Fig. [3] clearly demonstrate that £ is a lin- 
ear function of H 2 , without torsion, indicating a constant 
k = 0.69 eV. With the inclusion of the torsion term the 
resulting n increases with the radius of the nanotube. For 
large nanotubes, with small torsion angles, the value of 
0.69 eV found from the phonons is recovered. The value 



of 0.82 eV from Ref. [M] was indeed determined from a 
nanotube with a radius of approximately 11 A. 




Figure 3: Elastic energy per unit area, £, versus the curvature 
squared, H 2 , for nanotubes of varying radius. The slope of this 
curve determines the bending rigidity re. The inset shows the bend- 
ing rigidity determined for individual nanotubes. 



5. Nanotubes and multilayer graphene 

For completeness, we show in Fig. [4] also the phonons 
of a (10,10) nanotube that can be compared to Refs. [71 
|3"2"] . For nanotubes the differences between models are 
enhanced due the complex folding of the bands. 

Lastly, we examine the phonons of n-layer, AB stacked, 
graphene, going from a single graphene layer towards 
graphite. In this process, the low-frequency ZA mode 
splits into n optical sub-branches as shown in Fig. [5] 
while for all other branches the splitting is much smaller 
(^2 cm -1 ). As shown in the left panel of Fig. [5J 
the frequency of these 'breathing' modes at T for n- 
layer graphene are related to the longitudinal phonons of 
graphite along the T — A line at wavevectors 

_ 2irm . . . . 

<?"= , with (m = 0,...,n-l), (5) 

nc 

as if the modes were confined to an effective thickness of 
nc/2. This length is an interplanar distance larger than 
the actual thickness of the n-layer graphene. Interestingly, 
by extrapolating to a single layer, n = 1, we get an effective 
thickness equal to an interplanar distance, as suggested in 
Ref. 48J. Since the number and frequency of these low 
lying ZA modes is univocally determined by the number 
of graphene layers, their observation can be used for the 
characterization of multilayer graphene as a complement 
to the analysis relying on the 2D band done in Refs. [T^l 

EI]. 
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0.2 0.4 0.6 0.8 1 0.1 0.2 0.3 0.4 0.5 0.6 
Wavevector (jc/|T|) 

Figure 4: (a) Phonon dispersion from LCBOPII for a (10,10) car- 
bon nanotube; (b) The low frequency part of the phonon dispersion. 
|T| = 2.46 A is the lattice parameter along the nanotube axis. 



6. Conclusions 

Empirical potentials are desirable for their simplicity 
and transferability to calculate the phonon frequencies of 
complex systems. The phonon dispersions of graphite and 
graphene are an important test for the accuracy of these 
EIPs. We have shown that LCBOPII gives good results 
for graphitic crystals particularly in comparison to other 
EIPs. We have analyzed the reasons for the remaining 
discrepancies, suggesting that the potential could be im- 
proved considerably by modification of the long-range in- 
teractions. The quadratic ZA bending mode plays a key 
role in the graphene structure at finite temperatures and 
we have discussed how this fact might influence the fitting 
of force constants models to experimental values measured 
at room temperature. Lastly, we point out that multilayer 
graphene is characterized by several low frequency breath- 
ing modes at T, that are univocally related to the number 
of layers and could be used for their characterization. 
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